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Abstract 

Wave packets produced by localized disturbances play 
an important role in transition in three-dimensional 
boundary layers, such as that on a swept wing. Starting 
with the receptivity process, we show the effects of 
wave-space energy distribution on the development of 
packets and other three-dimensional disturbance 
patterns. Nonlinearity in the receptivity process is 
specifically addressed, including demonstration of an 
effect which can enhance receptivity of traveling 
crossflow disturbances. An efficient spatial numerical 
simulation method is demonstrated for the computation 
of these flows, allowing most of the simulations 
presented to be carried out on a workstation. 

Introduction 

Over the past ten years, technology for the 
prediction of transition to turbulence in boundary layers 
has progressed well beyond the use of quasi-parallel 
linear stability theory, with the recognition that much of 
the disagreement between theory and experiment were 
the result of physical effects beyond the scope of such 
methods. Among these are: receptivity - the process by 
which external disturbance energy is internalized into 
disturbances which grow in the boundary layer; effects 
of nonparallelism and surface curvature; nonlinear 
transfer of energy between disturbance modes; and 
effects of local inhomogeneities on the growth of 
disturbances. While numerical simulation technology 
has progressed to the point that full spatial DNS of the 
transition process is at least feasible and has been 
carried out for a few select cases, examination of the 
above effects across a broad parameter range is 
impractical due to the expense of such simulations. 

Many of these problems of interest may be cast 
in a linear or weakly -nonlinear framework, for which 
the modest spatial resolution requirements would render 
DNS a quite tractable tool for parametric studies. For 
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the most part, only the steady-state results of such 
simulations are of interest; the start-up transients, which 
are always present in simulations carried out in the 
time-dependent formulation, contain little useful 
information regarding the long-time asymptotic state of 
transition of a given flow. For instance, if a time- 
dependent simulation of a single-frequency wave is 
carried out from an impulsive start of the harmonic 
“driver”, then the transient will be observed as a 
leading-wave region which travels downstream; this 
region will consist of a broad spectrum of disturbances, 
initiated due to the step-function input. Each 
component wave will travel downstream at its 
characteristic group velocity, the slowest of which for 
the most part determines the length of the transient 
period. For many problems involving stationary or low- 
frequency disturbances, this period may be quite 
lengthy, and the simulation correspondingly expensive. 

However, if the assumption is made that the 
disturbance field is comprised of only a few select 
frequencies, then efficient use can be made of many of 
the algorithmic elements of full spatial DNS, while 
enabling the use of fast direct linear solvers. The 
advantage of such a methodology is that solutions may 
be obtained rapidly on workstation-level machines, 
enabling a researcher to perform the desired parameter 
studies of complex transition physics. This paper is an 
overview of the use of one such implementation over 
the past five years; studies carried out have included the 
generation of validation data for engineering transition 
prediction methods including receptivity-prediction 
methods, the detailed design of hybrid laminar flow 
control suction surfaces, studies of linear and weakly- 
nonlinear wave packets, and computations of the effects 
of streamwise and spanwise inhomogeneities on the 
evolution of disturbances. 
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Algorithmic Aspects 

What follows is a brief outline of the basic 
incompressible harmonic linear Navier-Stokes (HLNS) 
solver which was used in various forms to produce the 
results discussed below; a more complete description of 
the method will appear in a subsequent paper [1]. 

The first step in the formulation is expansion 
about a steady base flow: 

u = U b + eu’ 

where e is assumed (for now) to be a small parameter. 
After substitution into the three dimensional 
incompressible Navier-Stokes equations, equations of 
like order are separated. The 0(1) equations state that 
the base flow must satisfy the steady NS equations; 
however, we frequently use approximate solutions (such 
as solutions to the boundary -layer equations) for 
convenience, and for direct comparison with results 
from other methods. The 0(e 2 ) equations contain the 
nonlinear terms, and are used to compute weakly- 
nonlinear corrections to the HLNS solutions. The 0(e) 
equations, the linear disturbance equations, are of 
primary interest here. The use of a disturbance equation 
formulation is well-known to be required in 
computations of the initial stages of transition, as the 
accurate representation of the evolution of disturbances 
several orders of magnitude smaller than the base flow 
is essential. 

The following assumptions are then utilized in the linear 
disturbance equations: 

3uV3t = -ito u’ du73z = ip u’ 

the first of which reflects the single-harmonic 
assumption, the second being the assumption of 
spanwise (z) homogeneity. The latter is directly useful 
for oblique disturbances in a two-dimensional base 
state, and for disturbances evolving within the infinite- 
swept wing framework. For more complex flows, a 
Fourier integral method is used to represent the 
spanwise dependencies, as will be discussed later. 

The above process results in a complex two- 
dimensional equation set with two parameters: (O - 
related to the frequency of the disturbance, and P - the 
spanwise wavenumber. The streamwise (x) and normal 
(y) directions are discretized using relatively standard 
high-order methods: Chebyshev collocation is used in 
y, while fourth-order centered differences are used in x. 
These are natural choices given the solution scheme to 


be discussed below. Analytic mappings are used in both 
discretizations to improve resolution where required. 
Surface curvature is accounted for in the present method 
using the standard body-conforming thin-layer 
assumptions, which result in the simple addition of a 
few geometric factors to the Cartesian equation set; the 
Cartesian equations are smoothly recovered as the 
surface curvature becomes small. 

The resulting complex algebraic system of 
equations are in the form of a block pentadiagonal 
system; the blocks are of size 4N Y x4N Y , and are full 
due to the spectral discretization in the y-direction. The 
key to the usability of this method is the efficient direct 
solution scheme developed for the discretized equation 
set; the system is solved using standard recursion, but 
the recursion coefficients computed during the forward- 
sweep phase are stored on disk. Asynchronous I/O, in 
which the actual disk operations are carried out 
simultaneously which the computation of the next set of 
coefficients, speeds the overall throughput. In an 
average computation with N x = 1200 and N y = 51, a 
solution can be obtained in about 15 minutes on an SGI 
workstation; about 1.5 Gb of temporary disk storage is 
required. 

Boundary conditions are a sensitive aspect of 
this formulation; it was found during the early 
development of this method that some boundary 
conditions which are commonly used for time- 
dependent simulations can result in spectacularly poor 
and meaningless results when applied in the harmonic 
framework. This is due to two reasons: First, the 
condition number of the system is large to begin with, 
since the terms on the diagonal are relatively small. A 
boundary condition set which would have the effect of 
merely limiting the maximum time step allowable in a 
time-dependent simulation could render the direct 
solution so contaminated by roundoff error as to be 
useless. Second, the harmonic solution is long-time 
asymptotic; outflow boundary conditions (in particular) 
which build spurious error slowly over time result in 
solutions dominated by the effect of reflections when 
applied to this framework. The buffer-domain method 
[2], in which the momentum equations are smoothly 
parabolized to convect disturbances through the outflow 
boundary, was found to be robust and reliable here. 

Disturbance forcing is accomplished in many 
ways, depending on the particular physics being 
simulated. Regions or distributions of non-zero wall- 
normal velocity are used to simulate oscillating 
suction/blowing, for instance; formulations for 
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receptivity simulations result in either non-zero surface 
tangential velocity or field-forcing, as discussed below. 

Receptivity, the conversion of farfield 
disturbance energy into instability waves, results from 
the bilinear interaction between the unsteady signature 
of the farfield disturbance and a short-scale meanflow 
variation. In other words, a freestream disturbance of 
the proper frequency but long wavelength is scattered 
onto the short-wavelength instability spectrum by the 
local (steady) meanflow distortion of a surface 
inhomogeneity, such as a roughness element. See [3] 
for more details. The flowfield may be expanded as: 

U(x,y,t) = U c (X,y) + e„U„(x,y) + £ ac U ac (x,y)e + 
c.E-lUfx.yle”' 

where the coordinate M X” in U b (X,y) denotes slow 
variation in X. may be computed for very small 
roughness heights using the HLNS solver with co=0, and 
extrapolating the flow-iangency condition to y=0 to 
produce an inhomogeneous boundary condition on the 
tangential velocity: 

IL(X,0) = aU B (X,0) / Sy • h(x) 
where h(x) is the normalized roughness shape. 
Alternatively, the full distorted steady field 
corresponding to U b (X,y) + e w ll w (x,y) may be 
computed directly using a Newton-iterated nonlinear- 
solver variant of the HLNS method, where the 
roughness geometry is represented exactly using a 
shearing transformation. As will be described later, the 
local effect on the evolution of the recepted instability 
wave of this distorted base flow is the dominant 
nonlinear effect w hen the roughness height e w is finite. 

U.(x,y). the signature of the farfield 
disturbance, may be approximated in the case of 
acoustic receptivity by a Stokes wave; a more exact 
representation may be obtained through a solution of the 
linearized unsteady boundary- layer equations (LUBLE) 
[4]. Once U ac and U w are computed, their nonlinear 
interaction is used on the RHS of the HLNS solver for 
solution of the U wac problem. 

Since each (u),P) component is computed 
independently in the linear formulation, wave packets 
may also be easily constructed. For instance, the 
disturbance produced by an isolated roughness which is 
compact in the z-direction has energy distributed in the 
[3-plane; the Fourier integral of the bump shape provides 
the necessary wall -forcing as a function of p: 

m = ^]h(z)co s (/}z)dz 


The Fourier integral is discretized using a Chebyshev 
collocation quadrature formula: 

h(fi) = 'Z h (Zi) c °s(fiZj)Wi 

i = 0 

This discretization results in a set of { p, } at which the 
HLNS solutions are required. The packet is then 
constructed in the physical space via the inverse 
integral: 

U(x, y, z, t) = h(fi) 0 {x,y,/3,t) cos (/?z) dz 

« j^h{fi)V{x,y,fct)cos{/3,z) w, 

i = 0 

Results 

Crossflow Stability 

As mentioned in the Introduction, the effects of 
streamline curvature and nonparallelism on the 
development of disturbances in three-dimensional 
boundary layers has been a significant concern in the 
development and application of transition prediction 
methods, especially for swept wings. Extensive use of 
the subject HLNS method has been made in the 
evaluation of these effects and in the calibration of 
advanced transition prediction tools. 

To illustrate these effects, results are presented 
for crossflow disturbance growth in the boundary layer 
on the wing used in the ASU 45° swept-wing 
experiment [5], the pressure distribution for which is 
presented in Fig. 1 . The modest favorable pressure 
gradient over the first 70% of the wing gives rise to 
consistent growth of crossflow disturbances over a 
relatively narrow band of spanwise wavenumbers, and 
napthaline surface flow visualization indicated the 
presence of strong stationary (to = 0) crossflow 
disturbances. The effect of surface curvature on the 
evolution of one representative span wise- wavenumber 
disturbance in terms of the local growth rate evaluated 
using the maximum streamwise disturbance velocity is 
shown in Fig. 2; for comparison, the result from a 
standard quasi-parallel linear stability theory (LST) 
calculation is also shown. Note that the effect of 
surface curvature is to significantly reduce the local 
growth rate in this region near the leading edge; farther 
downstream, the two HLNS results merge as the 
curvature decreases. The large increase in growth rate 
shown by the HLNS results near the leading edge stem 
from the strong nonparallel effect, as the boundary layer 
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grows rapidly in the first 5%-chord arclength from the 
attachment line. The fact that the LST and HLNS 
results appear to agree after about x/c = .07 is 
completely fortuitous; results for other flowfields show 
significant disagreement between the very approximate 
LST and the virtually exact (for linear disturbance 
evolution) HLNS method. An additional result for this 
flowfield is shown in Fig. 3; in this figure, a comparison 


ASU Wing at a= —5° 



x/ c 

Figure 1 Pressure distribution on ASU swept wing 


between HLNS, PSE and LST results is made. 
Although the PSE results show some effect due to 
boundary-layer nonparallelism, this effect is still 
underestimated. Also shown in the figure is the effect 
of changing the streamwise placement of the 
suction/blowing strip which was used in the HLNS 
calculations to initiate the crossflow disturbance. 


ASU wing, /? = 900, co = 0 



X 


Figure 2 Stationary crossflow growth rate 


ASU wing, /?= 900 
w/ curvature 



Figure 3 Stationary crossflow growth rate, comparison of LST, PSE, and HLNS 
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Crossflow Receptivity - Packets 

Significant use of the HLNS method to predict 
the receptivity and evolution of crossflow disturbances 
in the presence of localized suction was made during the 
design and execution of a NASA / Boeing hybrid 
laminar flow control (HLFC) swept-wing experiment 
[6]; results from this test are proprietary, but overall 
indicated that standard LST transition prediction 
methods are quite unsatisfactory for use in swept-wing 
flows. 

The computation of the development and 
evolution of a packet of crossflow disturbances behind a 
roughness element near the leading edge of a swept 
wing will serve as an example of the use of the 
spanwise Fourier- integral method described above. A 
surface bump, circular in plan and Gaussian in cross- 
section, was represented using linear boundary 
conditions near the attachment line in a swept-wing 
boundary layer. Fig. 4 shows a plan view of isosurfaces 
of the streamwise component of disturbance velocity, 
red denoting positive and blue negative, with the lateral 
boundaries of the displayed computational region being 
parallel to the freestream direction and the inflow 
boundary (left) parallel to the leading edge of the wing. 
A particular pressure distribution was chosen for the 


infinite-swept mean flow which results in a rather wide- 
band growth of crossflow disturbances, and the small- 
scale roughness element which induces the disturbance 
packet in the computation creates initial energy in a 
range of scales which includes this growth band. The 
small variation of wave angle and growth direction 
across this band of growing crossflow disturbances is 
visible in the downstream evolution of the packet; 
disturbances with small spanwise wavelength appear on 
the downstream side of the packet, whereas those of 
longer wavelength propagate on the windward side. 

This effect is more apparent in Fig. 5, in which 
disturbance streamwise velocity from cuts at various 
chordwise locations is shown. 

Similarly, the initiation of stationary crossflow 
disturbances by perforated suction surfaces may be 
computed. In Fig 6 are shown isosurfaces of positive 
and negative disturbance velocity for flow of a swept- 
wing boundary layer over a spanwise-periodic array of 
suction holes; the 7 rows of holes in the simulation are 
evident in the figure. For this case, the spacings and 
angles of the rows of holes were chosen to produce near 
alignment of the holes along constant disturbance-phase 
lines, resulting in constructive interference of the waves 
from one row to the next. Simulations of this type have 
been used to study the performance sensitivity of 
perforate patterns for use in HLFC applications. 



Figure 4 Isosurfaces of disturbance velocity behind roughness element on swept wing. 
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When a stationary crossflow disturbance of a 
particular spanwise wavenumber p w evolves in a 
meanflow field that contains a small-amplitude slow 
spanwise variation, the disturbance is distorted. This 
distortion may be approximated, if the spectrum of 
the spanwise variation of the meanflow is confined to 
small wavenumbers, by the bilinear scattering of the 
P w into sidebands. For instance, if the meanflow 
variation possesses only a single spanwise component 
p o , then this bilinear interaction forces sideband 
disturbances of spanwise wavenumber p w ± P o . If P o 
is sufficiently small, these sideband disturbances will 
grow on their own, since their spanwise wavenumbers 
will lie in the same unstable range as the original 
disturbance. Thus, only a short stream wise region of 
meanflow variation is necessary to produce a 
disturbance field with rich spanwise content. 


10 


3 



Figure 7 Disturbance amplitude, effect of spanwise 
meanflow distortion on stationary crossflow 
disturbance 


The simple bilinear interaction is 
demonstrated by the example of a single crossflow 
disturbance, initiated by a roughness strip, evolving 
into an unstable mean boundary layer on which a 1 %- 
amplitude distortion with small spanwise 
wavenumber is imposed. In Fig. 7 is plotted the 
growth of the streamwise component of disturbance 
velocity; the initiation of the original wave (p w =4500) 
is apparent near the leading edge. The spanwise 
distortion of the meanflow occurs over the region .01 
< s/c; note that the scattered sideband disturbances 
grow rapidly to nearly comparable amplitude of the 
original disturbance. However, since the continued 
forcing of these waves is at a slightly different 
streamwise wavenumber than the disturbances which 
grow at the sideband spanwise wavenumbers (p = 
4186 & 4814), constructive and destructive 
interference occurs in the evolution of these waves. 
The reconstructed disturbance flowfield shows the 
expected spanwise variation in amplitude, as seen in 
Fig. 8. 



Figure 8 Contours of disturbance velocity, 
corresponding to Fig. 7 
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Nonlinear Receptivity 

The HLNS method has also been used 
extensively in the calibration of receptivity-prediction 
methods and theories, to be used in engineering 
amplitude-prediction tools for transition prediction. As 
described above, both linear and weakly nonlinear 
acoustic receptivity predictions are possible using 
HLNS methods. Two examples of receptivity of two- 
dimensional TS waves are shown here, the first 
resulting from a roughness strip, the second from a 
suction slot. The first case, shown in Fig. 9 in terms of 
receptivity coefficient against roughness height, 


Roughness Receptivity 



h (/xm) 

Figure 9 Roughness receptivity coefficient 


Distributed Receptivity 

Since the surface roughness distribution may 
be specified arbitrarily in the HLNS predictions, 
distributed roughness receptivity effects are easily 
simulated. Shown in Fig. 1 1 is the disturbance 
amplitude nriax y {U wac (x,y)} (without roughness height 
and acoustic amplitude parameters £ w and £„), 
corresponding to the experiment of Wiegel and Wlezien 
[9] in which a series of tape strips were placed at 
spacings corresponding to the wavelength of the TS 
wave excited at the acoustic-driver frequency in a flat- 
plate boundary layer. The receptivity saturation 
observed in that experiment, and predicted by an 
extension of the local receptivity theory [10], is 
apparent. 


indicates that the finite-Reynolds number local method 
of Choudhari, et al [7], is in complete agreement with 
the purely linear simulation result; however, significant 
deviation is seen between these results and the 
experimental results of Saric, et al [8]. The weakly 
nonlinear results, in which the effect of the roughness 
distortion on the local stability of the disturbance is 
taken account, shows much better agreement. 

Similarly, the local and HLNS results agree well for the 
suction-slot case shown in Fig. 10; however, the local- 
stability modification effect is even stronger for this 
case, and the weakly-nonlinear results deviate from the 
linear results significantly even for small suction rates 

Suction— Strip Receptivity 



Figure 10 Suction receptivity coefficient 


When the relative phases of waves induced by 
individual roughness elements in a distributed- 
receptivity situation do not precisely correspond, the 
receptivity is “detuned” and significantly weaker. Of 
course, for true distributed roughness on an aircraft 
surface, detuning would always be present and must be 
predicted, at least in the mean. An example of this 
detuning is shown in Fig. 12, in which the TS 
disturbance amplitude for acoustic receptivity at a 
particular frequency over a wavy wall is shown, first for 
the case in which the acoustic frequency corresponds to 
the wavelength of the wall, and second for the case in 
which the wavelength of the wall is 10% shorter. 
Evident is the destructive interference between the 
waves produced in the leading part of the wavy-wall 
region (denoted by the filled triangles) with the waves 
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produced later in the region. Work in this area of 
distributed-roughness receptivity, especially for the case 
of stationary crossflow receptivity, is continuing using 



0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 x!0 3 

K 

Figure 1 1 Disturbance amplitude, distributed 
roughness receptivity 


statistical models of roughness measured directly from 
surfaces of various levels of finish. 



0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 xlO 3 

K 

Figure 12 Disturbance amplitude, tuned and 
detuned wavy-wall receptivity 


Nonlinear Crossflow Receptivity 


In the case of traveling (i.e. non-zero 
frequency) crossflow disturbances, it was noted that the 
full distorted steady field corresponding to U b (X,y) + 
Ew U w (x,y) in the above formulation should also include 
the finite- amplitude stationary crossflow disturbance 
recepted by the roughness, as well as the nonlinear local 
distortion produced directly by the roughness. Since the 
streamwise wavenumber of the stationary disturbance is 
slightly different than that of the traveling disturbance, 
the interaction of the stationary disturbance with the 
Stokes' wave produces an additional receptivity 
mechanism that is detuned. This detuning results in an 
oscillatory amplitude increment over what would result 
from the interaction of the Stokes’ wave with the local 
roughness flowfield distortion. An example of this 
additional receptivity is shown in Fig. 13, for the case of 
a roughness strip with a particular spanwise wavelength, 
positioned near the leading edge of the ASU swept 
wing. Note that there are regions in which this 
additional receptivity mechanism adds an order of 
magnitude to the amplitude of the traveling crossflow 
disturbance, and that this increment is not a function of 
either £ w or 



Figure 13 Disturbance amplitude, nonlinear 
traveling crossflow receptivity 
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Harmonic Point Source 

In 1995, Watmuff [1 1 ] published results from 
very high resolution hot-wire measurements of the wave 
packet produced by oscillatory suction and blowing 
through a small hole in the surface under a flat-plate 
boundary layer. Ostensibly for studying the effect of 
discrete suction used for laminar flow control, the data 
provides an excellent test for approximate simulation 
methods. Initial comparisons with wave packets 
generated by summation of modes from linear stability 
theory, and with packets computed using the linear 
parabolized stability equations (PSE), showed that even 
for the small amplitude used in the experiment, some 
nonlinearity was present. This may be seen by 
comparing the isosurfaces of streamwise disturbance 
velocity from the experiment (Fig. 14) with those from 
the summation of 17 (symmetric) spanwise modes 
computed using the HLNS method (Fig. 15). Although 
the general wave pattern is captured by the linear 
computation, the experiment shows some detail on the 
centerline which is missing in the linear results; 
additionally, the experimental wave fronts are 
somewhat flatter than seen in the linear computation. 

Based on the conjecture that the disagreement 
was due to stabilizing meanflow distortion near the 
centerline of the packet, a nonlinear harmonic Navier- 
Stokes computation was carried out, in which only the 
primary, zero- frequency, and first harmonic temporal 
modes were included. The computation was done in a 
brute-force iterative manner, with the values of the 
nonlinear terms for each spanwise and temporal 
harmonic mode calculated from the previous global 
iteration and carried as forcing functions. These forcing 
functions were computed by exact decomposition for 
the temporal terms, and by reconstructing the total 
flowfield in successive streamwise planes, computing 
the nonlinear terms, and performing a spanwise FFT; 
the spanwise grid used for reconstruction used twice the 
number of spanwise modes, to dealias the computation. 
The results of this ongoing computation are shown in 
Fig. 16; the results near the centerline agree far better 
with the experiment, demonstrating the stabilization by 
the zero-frequency component; the computation does, 
however, suffer from lack of spanwise resolution. 

Conclusions 

A number of example applications of the 
HLNS method have been presented. It should be 
recalled that all of the results shown were computed on 
workstation- level hardware, and that broad parameter 


studies for many of the physical aspects of transition on 
swept wings discussed here have been carried out using 
a suite of codes developed around this method. 
Additionally, weakly nonlinear effects on disturbance 
growth, secondary instability, and disturbance scattering 
by meanflow distortions have also been computed using 
this method; these will be described in subsequent 
papers. 
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Figure 14 Isosurfaces of disturbance velocity, Watmuff HPS experiment 



Figure 15 Isosurfaces of disturbance velocity, linear HLNS computation 
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